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1. Introduction 

Perturbation theory offers a systematic way to improve theoretical predictions for 
any given infrared-safe observable. Depending on whether the expansion parameter 
in the series is ctg, or is as times a logarithm or a logarithm squared of some function 
of the kinematics of the hard process, one obtains a fixed-order result or a resummed 
result, respectively. Cross sections expanded up to a certain order in as and those 
resummed are relevant to complementary kinematic regions of the phase-space. It is 
therefore convenient to combine the features of these two expansions, by defining a 
matched cross section which is equal to the former or to the latter in the appropriate 
phase-space region. 

Fixed-order results are now generally available at the next-to-leading order (NLO), 
which corresponds to including in the cross section the coefficients of terms of order 
a\ (the leading order or Born level) and of order q;|+^. Thanks to the recent and 
rapid progress in the automated treatment of both the real and the virtual contribu- 
tions to NLO computations, it is realistic to assume that phenomenological results 
at this accuracy will become available in the next few years for all of the reactions 
of interest to the LHC physics programme. The situation is much less encouraging 
for fixed-order cross sections at (9(ag+^) (or NNLO), where only a handful of results 
are available, for very small final-state multiplicities. 

As far as resummed results are concerned, it is in general understood how to 
achieve a next-to-leading logarithmic (NLL) accuracy, which is equivalent to includ- 
ing terms proportional to a^ log*^" Q and to a" log*^"~^ Q, with k — 1,2 depending on 
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the nature of Q; some results are also known to next-to-next-to-leading logarithmic 
accuracy (NNLL). Unfortunately, resummed computations are in general technically 
complicated, observable-dependent, and error-prone; although for some observable 
classes a semi- automated algorithm (CAESAR [1]) is available, the overall situation 
is far less satisfactory than for NLO computations. For this reason, a very appealing 
alternative is that provided by Parton Shower Monte Carlos (PSMCs). Although 
formally PSMCs are equivalent to a LL-accurate resummation (which may become 
NLL in some corners of the phase space, for some PSMCs and subject to certain 
restrictions), in practice they are known to do much better than that, as compar- 
isons with data and with analytically-resummed results clearly show. Furthermore, 
PSMCs are fully flexible, give one the possibility of including hadronization mod- 
els in a consistent manner, and are the workhorses of experimental collaborations, 
thanks to their capabilities of simulating fully-realistic final-state configurations that 
can undergo detector simulations. 

As is well known, the approximations that form the core of PSMCs severely limit 
their predictive power in those phase-space regions (corresponding to multi-jet con- 
figurations) which are of interest for most of new-physics searches at colliders. These 
limitations can be alleviated by viewing PSMC predictions as resummed results, 
to be included with the NLO corrections to the relevant production processes into a 
matched cross section. The definition of a formalism for matching NLO computations 
and MC simulations has attracted a considerable amount of attention. There are now 
several proposals, but only two of them, namely MC@NLO [2] or POWHEG [3] , have 
made it to the stage of actually implementing several hadroproduction processes, and 
of being routinely used by experimental collaborations. 

The MC@NLO formalism requires the computation of the cross section predicted 
by the PSMC at C(as+^). Because of the structure of PSMCs, the non-trivial infor- 
mation of this computation is actually process-independent, and is contained in the 
definitions of the parton branchings; the process-dependent part is entirely factorized 
in the Born matrix element. Thus, one essentially has to perform one set of compu- 
tations (since typically initial- and final-state branchings are treated differently by 
the Monte Carlos) per PSMC, in order to be able to match an NLO computation 
with a parton shower simulation according to the MC@NLO approach. So far, these 
computations have been carried out for the case of HERWIG (see refs. [2, 4, 5]) and, 
more recently, for HERWIG++ [6]. In this paper, we present the first results relevant 
to the matching with PYTHIA 6.4. We limit ourselves to considering the case of 
initial-state branchings, and present some sample results for the Drell-Yan process. 

This paper is organized as follows. In sect. 2 we describe the various steps neces- 
sary for an MC@NLO matching with PYTHIA, from the definition of the underlying 
parton-level NLO cross section to the MC@NLO short-distance cross sections used 
for the generation of hard-subprocess events, to be showered by PYTHIA. We present 
sample results in sect. 3, and we give our conclusions in sect. 4. 
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2. MC@NLO 



2.1 NLO parton-level cross section 

The starting point for the construction of MC@NLO is that of writing the short- 
distance parton-level cross section according to the subtraction formalism of rcfs. [7, 
8] (which we shall call FKS subtraction henceforth). The default procedure in FKS is 
that of treating simultaneously (i.e. in one contribution to the partonic cross section) 
the two initial-state coUinear singularities, due to one given final-state parton being 
coUinear to the initial-state parton coming from the left or from the right. On the 
other hand, in FKS one can also treat these two singularities independently, by 
defining two separate contributions to the short-distance cross section, each of which 
corresponds to one of the initial-state singularities; this procedure is explained in 
detail in ref. [9]. As we shall discuss in the following, when an NLO computation 
is matched to PYTHIA according to the MC@NLO formalism, it turns out to be 
convenient to adopt the latter strategy (at variance with the case of HERWIG, where 
the simultaneous subtraction suffices). 

We shall assign the momenta entering the (real-emission) partonic subprocesses 
as follows: 



with V representing aW or a Z boson, and a, b, and c being QCD partons; Born-like 
processes are simply obtained from eq. (2.1) by removing c{k2) from the r.h.s.. Only 
soft and initial-state coUinear singularities are present in the processes of eq. (2.1). 
Hence, the S functions required for a separate treatment of the coUinear singularities 
are two and shall be denoted by 





S. 



(2.2) 



with the following properties: 




1, 
1, 



(2.3) 
(2.4) 



lim »S_ 



(2.5) 



lim S± 



^0. 



(2.6) 



The expectation value for any observable O will then be written as 



(O) = (O)^ + (O) 



(2.7) 
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where (see eqs. (4.14)-(4.16) of ref . [2]) 



= J dxidx2d(j)2 

ab 



0(2)5±(2) 



d(f)2 



(b) 



+ 0(1) ^ "^^"^ 



(c±) 



-0(1) 



1 



(c±) 
ab 



+ 



dE 



(sv) ' 
ab 



d(f)i 



Xi d4>idx 



ct 



-{0(i)5±(i),o(i)} 



(/) 

ab 



(2.8) 



As can be seen in ref. [2], the E terms in equation above are defined as the partonic 
short- distance cross sections, times the luminosity factors. 

In writing eq. (2.8), we have spUt (using the <S functions) the Born (E^^^) and 
soft- virtual (E^'^'''*) contributions into two terms, and have associated them with the 
corresponding real-emission contributions (see refs. [5, 9]). We have also used the 
coUinear hmits of eqs. (2.4) and (2.5), that imply that the S functions associated 
with the collinear counterterms are trivial. Finally, note that the coUinear remainders 
^ofe ii6ed not be multiplied by the S functions. 

Equation (2.8) must be further manipulated in order to use it in MC@NLO; in 
particular, one has to apply the so-called event projection which allows one to define 
a unique kinematic configuration associated with all counterevents (for a given real- 
emission configuration) - see sect. A.4 of ref. [2]. As discussed in that paper, although 
event projection is largely arbitrary, in the context of MC@NLO it turns out to be 
convenient to derive it from the behaviour of the Monte Carlo one interfaces to. We 
shall therefore discuss in the next subsection the behaviour of PYTHIA relevant to 
this issue. 

2.2 Event projection with PYTHIA 

The way in which PYTHIA deals with the simulation of V production is the follow- 
ing [10, 11]. First, the hard process is generated. This implies the generation of 
the two Bjorken x's entering such a process, which we shall denote by Ci and ^2 for 
the parton coming from the left and from the right respectively. The partonic cm. 
energy squared is 

So = (1(2-5, (2.9) 

with S the collider energy squared. In the case of single-V" production, sq = rriy, 
with mv the mass of V (or its virtuality if lepton pair production is considered), but 
eq. (2.9) is obviously valid regardless of the production process. PYTHIA then begins 
the showers. Choosing on statistical basis which leg emits "first" (such an emission 
is the only one that matters as far as MC@NLO is concerned), the shower variables 
are related to the momenta given in eq. (2.1) as follows: 

So 



(P1+P2)' 



(2.10) 
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(2.11) 



with i = 1, 2 for the emissions from the parton coming from the left and from the 
right respectively. The branching generated with eqs. (2.10) and (2.11) is such that 
the Bjorken x's associated with the initial-state partons of eq. (2.1) (which we shall 
denote by zi and Z2) are 

^l=Cl/^, Z2 = C2 (2.12) 

for an emission from leg 1, and 

zi = Cu Z2 = C2/Z (2.13) 

for an emission from leg 2. In the branching, the cm. energy of the Born-level 
subprocess is conserved (which, in the present case, is equivalent to the requirement 
that the virtuality of be a constant), and thus 

s = {pi+P2f = Z1Z2S = so/z. (2.14) 

An easy way to achieve this is that of imposing that the Bjorken x's of the Born-level 
subprocesses, Q, be separately conserved. 




Figure 1: Assignments of momentum fractions in the case of initial-state branchings as 
done in PYTHIA. Branchings from legs 1 and 2 are denoted in MC@NLO by + and — 
respectively. The blob represents the final state at the Born level. 



As shown in ref. [2], the event-projection procedure can be formally determined 
by constructing two "observables" that are conserved in the branching process. Ac- 
cording to what has been discussed above, a possible choice is: 

01 = CiC2^, (2.15) 

02 = ^log^. (2.16) 

It is worth noting that while eq. (2.15) is the same as in HERWIG, eq. (2.16) does 
not coincide with either of the HERWIG choices considered in ref. [2] ; the implication 
of this fact is that we obtain two different event projections, depending on whether 
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it is leg 1 or 2 that emits. We stress that eqs. (2.15) and (2.16), and the resulting 
manipulations we are now going to describe, do not depend on the nature of the 
Born-level final state (a V boson rather than - say - a three-jet configuration), 
owing to the fact that PYTHIA adopts the so-called s-approach [10] when doing 
initial-state branchings for all hadroproduction processes. In the context of the s- 
approach, the invariant mass (a/sq) of the final-state system at the Born level is kept 
constant during the branching. This allows one to formally replace V with the set 
of final-state particles at the Born level, and ki with the sum of their four momenta, 
which does not entail any changes to eqs. (2.9)-(2.14), from which we ultimately 
derive eqs. (2.15) and (2.16). The procedure described here is thus fully general. 




Figure 2: Assignments of momentum fractions according to the event-projection proce- 
dure, for the soft configurations. See the text for details. 

and is not restricted to \^-boson production. We recall that event projection is a 
way to re-write the soft and collinear counterterms in a pure-NLO computation, 
so as all counterterms associated with a given S contribution at the real-emission 
level (in the present case, or 5^rfE(^)) have the same kinematics. In order 

to achieve this, it turns out to be convenient to define the counterterm kinematics 
starting from a fixed real-emission kinematics. This situation is depicted in fig. 2 
for soft counterterms (the case of collinear counterterms is essentially identical, the 
only differences being in the assignments of the momentum fractions of the incoming 
partons, which we shall specify in what follows). Since the event projections we define 
here are ultimately motivated by the branching strategy of PYTHIA, the procedure 
of fig. 2 is by construction the inverse of that depicted in fig. 1. We shall return to 
this point in sect. 2.4. 

In order to actually obtain event projections, we use the master equations (A. 36) 
and (A. 3 7) of ref. [2]. For an emission from leg 1, these equations read in our case 

Oi(2(zi,Z2,02)) = Oi(l,(4t4t0?)), (2.1T) 
02(2(zi, Z2, 02)) = 02(l.(a;St xg, 4'^)) , (2.18) 
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and their analogues for the colhnear case (where one formally replaces s with c+ in 
the r.h.s. of eqs. (2.17) and (2.18)). Note that, having treated separately the two 
colhnear singularities in the NLO cross section, the case c— need not be considered 
when studying branching from leg 1). We find: 



X 



X 



(s) _ 
1+ " 
(c+) 
1+ 



Z Zi 



_ («) / 

— Xi^/Z , 



X. 



X 



(^) _ 

2+ - 

(c+) 
2+ 



Z2 



— X. 



is) 
2+ ■ 



For an emission from leg 2, we find instead 



x\i = 

(C-) 



X 



X 



1- ) 



(c-) 



ZZ2 , 
(s) 



X 



Iz. 



(2.19) 
(2.20) 

(2.21) 
(2.22) 



At this point, following what was done in sect. 4.4 of ref. [2], one uses eqs. (2.19) 
and (2.20) to perform the event projection transformation of (O)^, and eqs. (2.21) 
and (2.22) to deal with {0)_. Once event projection has been achieved, the MC@NLO 
short distance cross section can be defined by including the MC subtraction terms. 
We shall do this in the next subsection. 

2.3 MC subtraction terms 

At the NLO in as, the cross section resulting from PYTHIA is 



da 



MC 



+ da^-^ 



MC 



MC 



(2.23) 



where the two terms on the r.h.s. account for emissions from leg 1 and 2 respectively. 
They read: 



MC 



ahc 



E dC.dC^^da^iiCiPu C2P2)f f Pi?(.)/r)(Ci/^)/f ^HC2), (2.24) 



E^CiciC2^ciaiJ(0Pi,C2P2)^^i^:)(.)/f^)(Ci)/f^HC2/.). (2.25) 



ahc 



'27r 



ti z 



Following what was done in sect. A. 5 of ref. [2], we can manipulate the equations 
above to render them suitable for an integration together with the NLO parton-level 
cross section. First of all, one notes that the following identifications are vahd owing 
to the construction of event projections: 



in eq. (2.24), and 



C,l _ X^j^ , (,2 — X2j^ , 



C,l — Xi_ , C,2 — 2^2- ) 



in eq. (2.25). We can therefore perform a change of integration variables: 



MC 



dzidz2dT, 



ab 



ab 



MC 



(2.26) 
(2.27) 

(2.28) 
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with 



d^ab 



(+) 



MC 



ah 



MC z d{zi, Z2) 



MC 



MC 



and 



da 



(+) 



^^e(.)At™(4'>:,4'-;i'2)e, 



MC 



da 



(-) 



^ mil 



is) 



Here, we have introduced the rescaled virtuahties 



U = -ti/m 



V 



(2.29) 
(2.30) 

(2.31) 
(2.32) 

(2.33) 



Furthermore, we have exphcitly included the factors ©±, which are related to the 
choice of the maximum virtuality allowed during shower evolution. We have e.g. 



\ti\ < my 



l-3x 
l-x 



(2.34) 



with y being the relevant angular variable associated with the FKS parton in the 
FKS subtraction (here the cosine of the angle between ^2 and pi for emissions from 
leg i), and 1—x the (normalized) energy of the FKS parton. In the case of PYTHIA, 
at variance with HERWIG, x = z. 

2.4 MC@NLO short-distance cross sections 

At this point, we have all the ingredients needed to generate hard-subprocess events 
according to the MC@NLO formalism, which will be subsequently showered by 
PYTHIA. Recalling that we denote by e''"^ the contributions S*^") to the NLO parton- 
level cross sections after event projection, event generation will be performed follow- 
ing the procedure outhned in sect. 4.5 of ref. [2] (see also sect. 3 of ref. [5]), using 
the following integrals: 



r(±) 



ab 

/ dzidz2 



dY!-^^ 

dZidZ2d(t)2[ S±—^ 



dT. 



■(±) 

ab 



(2.35) 



ab 



+ 



-5- 



=(6) j-^{sv) 



^ f d^ab , dT,, 
Ii \ d(f)i 



+ 



-'ab 



d(f)i 



d(t)2 



+ 



+ 



ct 



d(t)2 



ab 



Xj \ d(f)\dx 



ab 



d(f)idx 



ct 



(2.36) 
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We remind the reader that all integrals defined in these equations are separately 
finite, and that 



(^tot — It 



(+) 



j(+) , j(-) 



(-) 



(2.37) 



is an exact equation (with atot the fully-inclusive NLO rate). 

Although eqs. (2.35) and (2.36) solve the problem of the generation of hard- 
subprocess events, a further simplification is possible. We start by considering the 
integrals that are used to generate H events. These two integrals have to be 
computed separately. However, one observes that, for a given choice of (-^i, 2^2, ^2), 
the kinematics associated with these two contributions are actually identical, and 
the integrals can therefore be computed together. This is equivalent to using the 
following integral for the generation of H events: 



j(+) , j(-) 



=i:/-.-*(|fL-|:f^ 



(2.38) 



The first term in the integrand has been simplified thanks to eq. (2.3). For the com- 
putation of the second term in the integrand we stress that, as eq. (2.38) explicitly 
indicates, Zi are the integration variables. Therefore, when we evaluate eqs. (2.31) 
and (2.32), the variables x['^2 have to be computed according to eqs. (2.19) and (2.21). 
From the physical viewpoint, this implies that the two Born matrix elements ap- 
pearing implicitly in E^^^ are evaluated in two different kinematics configurations, 
that eventually result after the first branching in the same real-emission kinematics 
{zi, Z2, (t>2}- This situation is precisely the one depicted in fig. 2. 

A simphfication analogous to that of eq. (2.38) is also possible in the case of § 
events, although the argument is slightly more involved. If one fixes {zi^ Z2, 4>2) (i-e., 
the real-emission kinematics), the kinematic configurations associated with and 
/g~'* in eq. (2.36) are different, owing to the fact that the two integrals are computed 
using different event projections. Formally, one can introduce two mappings, Pj^^g 
and P]H[J.§, representing these event projections; for a given real-emission kinematic 
configuration 2, the two Born-hke configurations associated with 7^"*"^ and can 



be denoted by 



1 + 



2 



)(-) 



(2.39) 

respectively; 1+ and 1~ are the hard configurations to be showered by PYTHIA. 
The mappings are pictorially represented in fig. 2 by the two thick arrows, 

with 2, 1"*", and 1~ the configurations depicted in the upper part, lower left corner, 
and lower right corner of that figure respectively. On the other hand, one can fix a 
Born-like configuration 1 and, using the inverse of the maps "Pe^g (which we denote 
by "Ps^^jj), work out the real-emission configurations to be used in the computation 
of the integrand of /g^^ : 



)+ 



)(+) 



>{-) 



(2.40) 
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As we have already stressed, although there is an ample freedom in choosing "P^jj^g, 
it is best to adopt a form motivated by what the PSMC does when branching. In 
practice, it is therefore convenient first to obtain "Pg^^jj from the PSMC, and then (by 
inverting them) V^^. In the case of PYTHIA, Pg^H amount to performing a boost 
of the Born-level four momenta in the transverse direction, to balance the transverse 
momentum of the parton produced in the branching (c(/c2) in eq. (2.1)), followed by 
a boost in the longitudinal direction, to e.g. the lab frame. The information on the 
latter is equivalent to the assignments of the momentum fractions of the incoming 
particles. The prescription for the transverse boost given above is obviously trivial 
in the case of V production (since ki^ = — A;2t), but can be applied to final states 
with arbitrary multiplicity. Thus, what done here is also valid for final states more 
complicated than single V. 

The two procedures related to eq. (2.39) and (2.40) are equivalent, since the inte- 
grals are performed over the whole phase space (possibly subject to kinematic restric- 
tions, which are however identical in the two cases). Equation (2.40) corresponds to 
the situation depicted in fig. 1, except for the seemingly different assignments of the 
momentum fractions of the incoming particles. These differences are however imma- 
terial, since such fractions are simply integration variables which can be manipulated 
and renamed at will. We shall now proceed to perform such manipulations, starting 
from a change of variables in eq. (2.36) 



(^1,^2,02) 



(2.41) 



Furthermore, consistently with eq. (2.40), the integrands have to be computed with 
the suitable kinematic configurations, 2+ or 2~. We can take these two operations 
into account by the formal replacements 



2± 



d{zi,Z2) ^(a), 



dixit xii) 



-"ab 



2(^^.(42))) • 



(2.42) 



It is easy to realize that, in the case of the Born, soft- virtual, and soft counterterms 
contributions, we have by construction 



-'ab 



(a) 
ab 



(2.43) 



since for these terms the jacobian appearing in eq. (2.42) exactly compensates the 
one in the definition of the S^^^ terms, see eq. (4.18) of [2]. In general, this is not true 
for contributions with a purely coUinear structure (such as the coUinear counterterms 
to the real-emission matrix elements). Equations (2.41), (2.42), and (2.43) allow us 
to rewrite 



r(±) 



E 

ab 



- S- 



+ 



+ 



dE 



ab 



d(f)i 



+ 



#2 
1 



+ 



ct C?02 
(c±) 



MC 



dE 



ab 



Xj \d(f)idx 



(c±) 
ab 

d(f)idx 



ct 



(2.44) 
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At this point, the two integrals in eq. (2.44) correspond by construction to the same 
Born-hke kinematic configuration (to be showered by PYTHIA), and can therefore 
be integrated together. It is thus convenient to rename the integration variables 



{x 



is) 



X 



2±> ) 



(Cl,C2,02) , 



and use 



r(+) 



(-) 



(2.45) 



(2.46) 



for the generation of S events. One observes that in the integrand of eq. (2.46) there 
will be a term 



d4> 



1 

2i 



dcf) 



d4> 



(2.47) 



where we have made use of eq. (2.3). The same kind of simplification will occur in the 
sum of the two soft counterterms. As far as the coUinear counterterms are concerned, 
we stress that the S± terms that multiply them are equal to one (by construction of 
the S functions - see eqs. (2.4) and (2.5)). We arrive therefore at the following form: 



J2 / dCidC2dh 



ab 

1 



dT. 



if) 

ab 



dT. 



(b) 
ab 



+ 



dT. 



#2 

(sv) 



ab 



d(f)i 



45: 



L=± 









\d(t)idx 


d(pidx 


ct^- 



(2.48) 



Equations (2.38) and (2.48) are our final expressions, used for the generation of H and 

S events respectively. The net result of the various manipulations carried out in this 
section is that it is still possible to treat simultaneously the two initial-state collinear 
singularities, as was the case when matching with HERWIG. We point out that the 
fact that eqs. (2.38) and (2.48) do not depend upon the S functions is not a property 
of V-boson production, but applies as well to all processes whose Born-level final-state 
particles are all colour singlets^. In fact, for this independence of <S functions to occur, 
we only need eq. (2.3) to hold, which is true in the cases just mentioned (see ref. [9] 
for further details on the construction of S functions). Clearly, in order to make 
use of eq. (2.3), iS+ and <S_ need be computed for the same kinematic configuration. 
This is in fact what happens, in spite of the fact that in the intermediate steps of the 
event-projection procedure we had to treat separately the + and — contributions. For 
EI events this is trivially true by construction - both I^'' in eq. (2.35) arc associated 
with the kinematic configuration denoted by 2 in eq. (2.39) (i.e. the upper part of 
fig. 2). On the other hand, the short-distance cross sections for S events we started 



^It is also possible, but not mandatory, to treat in this way processes in which all strongly- 
interacting final-state particles are massive, as e.g. in tt production. See ref. [4] for an explicit 
example. 
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from, eq. (2.36), are associated with two different kinematic configurations, denoted 
by 1^ in eq. (2.39) (i.e. the lower part of fig. 2). Thanks to the procedure described 
above, we have manipulated eq. (2.36) precisely to be able to associate the two ± 
contributions with the same kinematic configuration, denoted by 1 in eq. (2.40) (i.e. 
the upper part of fig. 1). Owing to the properties of event projections that we have 
discussed above, the whole procedure can obviously be carried out for any final-state 
multiplicity, with 2 and 1 formally replaced by the real-emission and Born-level 
configurations respectively. Finally, we stress that the manipulations performed here 
will also be valid in the case of strongly-interacting particles in the final state. The 
only difference is that, in such a case, eq. (2.3) will not hold any longer, and therefore 
the quantity iS+ + S- will appear in eqs. (2.38) and (2.48) as a factor multiplying the 
real-emission, soft counterterms, soft-virtual, and Born contributions (on the other 
hand, the purely coUinear terms will be unchanged, owing to the fact that eqs. (2.4) 
and (2.5) hold regardless of the nature of the final state). Clearly, the contributions 
due to branchings of Born-level final-state strongly-interacting particles will be given 
by short-distance cross sections analogous to those of eqs. (2.38) and (2.48) (see e.g. 
sect. 3 of ref. [5]). The computation of such contributions is beyond the scope of the 
present work, and we postpone it to a future publication. 

We conclude this section by stressing that, due to the form of the shower variables 
used by PYTHIA, the soft limit of the sum of the MC subtraction terms coincides 
with that of the real-emission matrix elements, which was not the case for HERWIG. 
This implies that here we can set ^ = 1, where Q is the function introduced in 
sect. A. 5 of ref. [2] - we refer the reader to that paper for a discussion on this issue. 

3. Results 

In this section, we present sample results relevant to production in pp collisions 
at y/S = U TeV. Our aim is not that of performing a phenomenological study, but 
rather that of presenting a few control plots that show that the matching of the NLO 
results for single vector boson production with PYTHIA according to the MC@NLO 
formalism works as we expect. This is non trivial, given the differences between 
HERWIG and PYTHIA, which in turn result in different short-distance MC@NLO 
cross sections, as discussed previously. We set mw — 80.4 GeV, = 2.14 GeV, 
and use CTEQ6.6 [12] PDFs. Our default scale choices arc fip = fif^ = rrij,, where 
ttIt is the transverse mass of the W. When reconstructing jets, we adopt the /cT-jet- 
finding algorithm of ref. [13], with y^ut = (10 GeV)^. 

In each of the plots we present in figs. 3-5, we display three histograms. The 
solid (black) histograms are the results of MC@NLO with PYTHIA (which we shall 
call MC@NLO/PY) i.e. what has been computed in this paper. The dotted (red) his- 
tograms are the results of MC@NLO with HERWIG (which we shall call MC@NLO/HW) 
Finally, the dashed (blue) histograms are the results obtained with PYTHIA stan- 



12 



dalonc. When using PYTHIA for showering hard events, we set MSTP(81)=0 and 
MSTP(91)=0, which corresponds to switching off multiple interactions and primor- 
dial kj, respectively. We switch off matrix element corrections by setting MSTP(68)=0, 
which also forces the maximum virtuality in the shower to be equal to the vector 
boson mass, when standalone generation is performed. On the other hand, hard 
subprocess events generated by MC@NLO/PY are given to PYTHIA in the standard 
Les Houches format [14]; we set PARP(67)=1 and SCALUP=mvK to have the same 
maximum scale in the shower as in the standalone generation. In order to facili- 
tate the visual comparisons between MC@NLO/PY and PYTHIA, the results of the 
latter are rescalcd (by different factors, depending on the obscrvablcs considered). 
As has already been discussed at length in the literature, we expect MC@NLO to 
be identical in shape to the PSMC results in the regions of the phase space domi- 
nated by those large logarithms the PSMC is able to resum; this is the motivation 
for comparing MC@NLO/PY with PYTHIA standalone. Also, we expect MC@NLO 
to coincide, in shape and normalization, with the NLO results in the regions where 
hard emissions are dominant. Given that it has already been shown in the past 
that this is the case for MC@NLO/HW, here we compare MC@NLO/PY directly with 
MC@NLO/HW, rather than with the pure NLO results, since this will also give us the 
opportunity to observe the different behaviours of the two underlying PSMCs in the 
soft/collinear regions. As far as negative weights arc concerned, for W production 
at the LHC their fraction with MC@NLO/PY is about 0.6%, while in the case of 
MC@NLO/HW is about 8%. Since the phase-space parametrizations used in the two 
codes are identical, this large difference is essentially due to the different choices of 
shower variables made by the two PSMCs. 

In fig. 3 we consider the transverse momentum distributions of the boson 
(left pane), and of the hardest jet of the event (right pane). These distributions 
are interesting since in the large-p^ region they are dominated by hard emissions, 
whereas in the small-py one there are logarithmically-enhanced terms that need be 
resummed (the pure NLO results diverge there). The comparisons among the three 
results for the two distributions follow the same pattern, as we expect. At small prS, 
the shapes of MC@NLO/PY and of PYTHIA are identical. To show this clearly, the 
latter results have been rescaled so as their first bins coincide with those resulting 
from MC@NLO/PY. On the other hand, the MC@NLO/PY and MC@NLO/HW results 
are quite different in this pr region, owing to the different treatment of soft and 
collinear emissions by the two PSMCs. It is known that PYTHIA tends to give 
softer spectra than HERWIG, which is therefore what we expect to and do find 
with MC@NLO/PY and MC@NLO/HW. At large p^'s, MC@NLO/PY coincides with 
MC@NLO/HW - both in shape and in absolute normalization, which is again what 
we expect. 

In fig. 3 we also show the scale dependence of MC@NLO/PY, determined ac- 
cording to the following procedure. The renormalization-scale variations are de- 
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Figure 3: Solid (black): MC@NLO/PY. Dotted (red): MC@NLO/HW. Dashed (blue): 
PYTHIA standalone (rescaled), without matrix element corrections. Left pane: Pt of the 
W boson. Right pane: Pt of the hardest jet of the event. The insets show the fractional 
scale dependence of MC@NLO/PY, computed as described in the text. 



fined as the differences between the results obtained with (/i^j/i^,) = (771^,772^) (i.e. 
the default), and those obtained with {^ii,fip) = {fmT,mT), where / = 1/2 and 
/ = 2; these differences are computed bin-by-bin for all observables studied. Like- 
wise, for the factorization-scale variations one considers (yU^,/!^) = {rrij,, frriT)- The 
renormalization- and factorization-scale variations of like sign are then summed in 
quadrature, and the results are then summed to (for positive variations) or sub- 
tracted from (for negative variations) the default cross section. The ratios of the two 
predictions obtained in this way over the default cross section are displayed in the 
insets of fig. 3 as the upper and lower bounds of the shaded areas, for the two observ- 
ables considered there. In spite of the lack of statistics in the high-p^ tails^ the trend 
is clear: the fractional scale dependence grows from about ±5% at low p^'s to about 
±10% at large p^'s. This has to be compared with the pure-NLO result (not shown 
here) that features a decrease in the scale dependence, from ±15% at low pr^s to 
±10% at large Pt's; this behaviour results from a dependence on fip almost constant 
w.r.t. Pt, and a dependence on fin decreasing with (due to the running of as). The 
scale dependences of MC@NLO/PY and of the pure-NLO predictions are therefore 
consistent at large p-^'s. This is what we expect, since there the MC@NLO/PY result 
is basically coincident with the NLO one, the effect of the shower being negligible on 
these inclusive variables. The size of the scale dependence is also compatible with 
the fact that, in fixed-order perturbation theory, 0{as) is actually the first order 

^MC@NLO outputs unweighted events. The plots presented here have been obtained with 5-10^- 
event samples, and therefore less than five hundred events have p^'s larger than a few hundreds 
GeV. 
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that contributes to Pt > 0. This observation therefore apphes also to low p^'s in 
the case of the pure-NLO predictions, but it does not in the case of MC@NLO/PY. 
In fact, MC@NLO/PY fills the low- and intermediate-pr regions mostly through the 
showering of § events. The shower, however, determines only the kinematics of the 
final-state configuration (i.e., the Pt of the W and of the hardest jet here), but the 
weights are given by the short- distance cross section of eq. (2.48), and therefore re- 
ceive both 0(0°) and 0{ag) contributions. This implies that in the \ow-pj, region 
one expects MC@NLO/PY to have a scale dependence smaller than that of the tree- 
level 0{as) term alone, which is what we observe. We conclude this discussion by 
noting that other definitions can be given of the uncertainties associated with mass 
scales (e.g., variations may be summed linearly rather than in quadrature), without 
this changing the pattern found here. We also point out that the scale dependences 
of the rapidity observables which we shall discuss below are rather featureless (i.e., 
scale variations are independent of rapidities), and therefore they will not be shown 
in what follows. 



1 I ' ' ' ' I ' ' ' ' I ' ' ' ' I ' ' ' ' I ' ' ' ' I ' ' ' ' 




Figure 4: Same as in fig. 3. Left pane: rapidity of tlie W boson. Right pane: rapidity of 
the hardest jet of the event. 

In fig. 4 we consider the rapidity distributions of the boson (left pane), 
and of the hardest jet of the event (right pane). In the absence of any cuts in 
Pt, the boson rapidity is an inclusive variable unaffected by large logarithms, and 
both NLO computations and MC simulations should predict it relatively well. We do 
indeed see an overall consistency among MC@NLO/PY, MC@NLO/HW, and PYTHIA 
standalone. Differences among the three predictions are larger in the case of the 
rapidity of the hardest jet, since this is a less inclusive variable w.r.t. the W rapidity, 
and is also more sensitive to hadronization corrections. The PYTHIA standalone 
results have been rescaled by the NLO K-factor, at^Lol^Lo^ before any cuts are 
applied and jets are reconstructed. 



15 




-6 -4 -2 2 4 6 -6 -4 -2 2 4 6 



y(w)-y(ji,rd) y(w)-y(jtrd) 

Figure 5: Same as in fig. 3, for the difference in rapidity between the W boson and the 
hardest jet of the event. Two cuts on the px of the jet are considered. 

Finally, in fig. 5 we consider the difference between the rapidities of the 
boson and of the hardest jet, for two different cuts on the transverse momentum of 
the hardest jet. The PYTHIA standalone results have been rescaled in the same way 
as in fig. 4. This observable and its analogues (with yw replaced hj ys, S being 
the system emerging from the hard process at the Born level - e.g. the ti pair in 
top-pair production) has attracted some attention in the past (see e.g. ref. [15] for 
a recent discussion), owing to the fact that MC@NLO/HW has a different behaviour 
around ys ~ yj — w.r.t. the underlying NLO computation - the former being 
flatter than the latter, or having a dip, depending on the nature of the system S. 
The "flatness" or the presence of a dip is a feature of MC simulations, and more 
specifically is a consequence of the choices made for initial conditions, as was done 
here for PYTHIA in eq. (2.34). It has to be stressed that this is true for both HERWIG 
and PYTHIA, as the dashed histogram on the right pane of fig. 5 shows, and therefore 
has nothing to do with the presence of dead zones^ as such in HERWIG. In PYTHIA, 
the choice of initial conditions is much more flexible than in HERWIG, and one can 
make the dip of fig. 5 disappear by choosing a large-enough scale as the maximum 
virtuality allowed for the shower. In doing so, however, one may extend the collinear 
approximation, that is the core of both HERWIG and PYTHIA, outside its proper 
range of validity, and we consider the choice made in eq. (2.34) a sensible one. EI 
events in MC@NLO, or matrix element corrections in HERWIG and PYTHIA, will 
give a substantial contribution to the region ys — yj ^ (see the solid and dotted 
histograms in fig. 5), where the predictions will become closer (w.r.t. those of MC 
simulations without matrix element corrections) to the pure-NLO results. Differences 
will in general remain, that can be formally ascribed to effects beyond NLO. 

^Negative weights in MC@NLO are also not an issue. This is a fortiori true in MC@NLO/PY, 
where they are basically negligible. 
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In order to further this argument, we consider an extension of the choice made 
in eq. (2.34). Namely, we parametrize the maximum virtuahty allowed in the shower 
in terms of a number / introduced as follows: 



\ti\ < ifmvY 



eh/ 



l-(l + 2/2)x 



X 



(3.1) 



The plots shown so far have therefore been obtained with / = 1. In fig. 6 we consider 
again the difference in rapidity of fig. 5, with / = 1/2 (dotted red histograms) 
and with / = 2 (dashed blue histograms), together with our default choice / = 1 
(solid black histograms); we point out that for the values of / considered here the 
fraction of negative weights is basically a constant. Two additional Pt cuts are also 
considered on top of those of fig. 5. The results obtained with PYTHIA feature a 
very large sensitivity to the choice of /, while those obtained with MC@NLO/PY are 
fairly stable. One may be tempted to use the results of fig. 6 in order to "tune" 
the parameter /, and obtain a PYTHIA prediction in decent agreement with that 
of MC@NLO/PY. However, this is an a posteriori procedure, which can be justified 
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Figure 6: Rapidity differences between the W and the hardest jet of the event, with 
PYTHIA (left pane) and MC@NLO/PY (right pane), for different values of / in eq. (3.1). 
Four cuts on the pr of the jet are considered. 



phenomenologically, but which implies that any possibility is given up of a sensible 
estimate of the theoretical uncertainties affecting the observable considered here; 
also, such a procedure is in general observable-dependent. Indeed, the left pane of 
fig. 6 shows that, if no information on 2 — )■ 2 matrix elements is included, this rapidity 
difference basically cannot be predicted hj Monte Carlos at moderate and large p^'s. 
This is not surprising, given that one is attempting to use the coUinear approximation 
outside its range of validity: the correct result can be recovered (which is equivalent 
to choosing the argument of a logarithm so as its numerical value coincides with a 
given constant), but only in a heuristic way. The situation improves if the correct 
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information on matrix elements is used, as in MC@NLO/PY. This is analogous to 
what happens when one varies the renormalization and factorization scales, at the 
LO and the NLO levels. As in the case of scale variations, extreme values for / 
will lead to problems; however, / = appears to be a safe choice, allowing one 
to realistically estimate theoretical uncertainties. We conclude by mentioning that 
other observables (such as those shown in fig. 3 and 4) display the same pattern of 
dependence on / as the observable in fig. 6. We shall discuss this issue further [16] by 
considering also the case of Higgs production, where the effects are more pronounced. 

4. Conclusions 

We have presented the construction of the matching between an NLO QCD compu- 
tation and the virtuahty-ordered PYTHIA Monte Carlo, according to the MC@NLO 
formahsm. We have hmited ourselves to considering only the case of initial-state 
radiation, and applied the formalism to the study of W hadroproduction. Owing to 
the different structures of HERWIG and PYTHIA, the short- distance MC@NLO cross 
sections used to generate hard events are different in the two cases. However, the 
FKS subtraction, which is the method used in MCCilNLO for dealing with infrared 
singularities, needs no modifications, and is able to treat both MCs. Likewise, no 
changes are needed in the MC@NLO formalism, and the matching with PYTHIA is 
achieved by performing a process-independent calculation. 

For the process considered in this paper, the fraction of negative weights in 
MC@NLO/PY is much smaller than that in MC@NLO/HW. Prom the physical view- 
point, the pattern of the comparison between MC@NLO, pure-NLO, and PSMC 
results is the same for MC@NLO/PY as for MC@NLO/HW — MC@NLO shows the 
same behaviour as the NLO or MC where either one is most reliable, with a smooth 
transition between the hard and soft-collinear emission regions. 

Although the process studied in this paper is particularly simple, the formulae 
given here will be basically sufficient for performing the matching in the case of more 
complicated reactions (and limited to initial-state branchings only). Therefore, the 
results presented here are the first step towards the construction of MC@NLO/PY 
for generic processes, where also final-state radiation is present, and towards the 
extension of this formalism to the pr-ordered version of PYTHIA. 

Acknowledgments 

We would hke to thank Torbjorn Sjostrand, Peter Skands, and Bryan Webber for 
helpful discussions and comments on the manuscript. The work of P. T. is supported 
by the Swiss National Science Foundation. 



18 



References 

[1] A. Banfi, G. P. Salam and G. Zanderighi, JHEP 0503 (2005) 073 [arXiv:hep- 
ph/0407286]. 

[2] S. Prixione and B. R. Webber, JHEP 0206 (2002) 029 [arXiv:hep-ph/0204244]. 
[3] P. Nason, JHEP 0411 (2004) 040 [arXiv:hep-ph/0409146]. 

[4] S. Prixione, P. Nason and B. R. Webber, JHEP 0308 (2003) 007 [arXiv:hep- 
ph/0305252]. 

[5] S. Prixione, E. Laenen, P. Motylinski and B. R. Webber, JHEP 0603 (2006) 092 
[arXiv:hep-ph/0512250]. 

[6] S. Prixione, S. Latunde-Dada, F. Stoeckli, P. Torrielli, in preparation. 

[7] S. Prixione, Z. Kunszt and A. Signer, Nucl. Phys. B 467 (1996) 399 [arXiv:hep- 

ph/9512328]. 

[8] S. Prixione, Nucl. Phys. B 507 (1997) 295 [arXiv:hep-ph/9706545]. 

[9] R. Prederix, S. Prixione, P. Maltoni and T. Stelzer, JHEP 0910 (2009) 003 
[arXiv:0908.4272 [hep-ph]]. 

[10] T. Sjostrand, S. Mrenna and P. Z. Skands, JHEP 0605 (2006) 026 [arXiv:hep- 
ph/0603175]. 

[11] G. Miu and T. Sjostrand, Phys. Lett. B 449 (1999) 313 [arXiv:hcp-ph/9812455]. 

[12] P. M. Nadolsky et al, Phys. Rev. D 78 (2008) 013004 [arXiv:0802.0007 [hep-ph]]. 

[13] S. Catani, Y. L. Dokshitzer, M. H. Seymour and B. R. Webber, Nucl. Phys. B 406 
(1993) 187. 

[14] E. Boos et al., arXiv:hep-ph/0109068. 

[15] S. AlioU, P. Nason, C. Oleari and E. Re, JHEP 0807 (2008) 060 [arXiv:0805.4802 
[hep-ph]]. 

[16] P. TorrieUi, in preparation. 



19 




1 2 3 



Logio(PTW/GeV) 




1 2 3 



Logio(PT(jhrd)/GeV) 



